Distinguishable and Exchangeable Dyads: Multilevel Modelling

Cross-Sectional and Intensive Longitudinal APIM and DIM

Pascal Küng

2025-10-20

Overview

  • Distinguishable Dyads
    • Cross-Sectional APIM
    • Longitudinal APIM
  • Exchangeable Dyads
    • Cross-Sectional APIM
    • Cross-Sectional DIM
    • Equivalence between APIM and DIM
    • Longitudinal DIM
    • Longitudinal APIM

Loading Libraries and setting CmdStan Backend

library(tidyverse)
library(brms)
library(bmlm)
library(easystats)
library(DiagrammeR)
library(glue)
library(DHARMa)
source(file.path('00_R_Functions', 'ReportModels.R'))
source(file.path('00_R_Functions', 'PrettyTables.R'))
source(file.path('00_R_Functions', 'PrepareData.R'))

set.seed(123)

options(
  brms.backend = 'cmdstan',
  brms.file_refit = 'on_change'
)

Distinguishable Dyads

Distinguishable Dyads - Cross-Sectional APIM

Distinguishable Dyads - Data: Simulated Dyads

df <- readRDS('Embed/simulatedCrossSectionalDyadicData.rds')
print_df(head(df))
userID coupleID gender communication satisfaction
1_1 1 1 4.851107 5.214440
1_2 1 2 6.004533 5.812518
2_1 2 1 5.881321 7.467339
2_2 2 2 6.433723 7.256714
3_1 3 1 4.283971 6.047375
3_2 3 2 5.516060 5.314505

Distinguishable Dyads - Cross-Sectional APIM

Distinguishable APIM

Distinguishable Dyads - Preparing Data

df_apim <- df %>%
  # Reshaping to Actor-Partner format (4-field)
  reshape_dyadic_data(
    person_id = 'userID',
    dyad_id = 'coupleID',
    vars_to_reshape = 'communication'
  )  %>%
  mutate(
    # Optional: grand-mean centering
    communication_actor_gmc = communication_actor - mean(communication_actor),
    communication_partner_gmc = communication_partner - mean(communication_partner),
    
    # Create Dummy-Variables for male and female
    is_male = ifelse(gender == 1, 1, 0),
    is_female = ifelse(gender == 2, 1, 0),
  ) 

print_df(head(df_apim))
userID coupleID gender is_male is_female communication satisfaction communication_actor communication_partner communication_actor_gmc communication_partner_gmc
1_1 1 1 1 0 4.851107 5.214440 4.851107 6.004533 -0.1553604 0.9980657
1_2 1 2 0 1 6.004533 5.812518 6.004533 4.851107 0.9980657 -0.1553604
2_1 2 1 1 0 5.881321 7.467339 5.881321 6.433723 0.8748537 1.4272563
2_2 2 2 0 1 6.433723 7.256714 6.433723 5.881321 1.4272563 0.8748537
3_1 3 1 1 0 4.283971 6.047375 4.283971 5.516060 -0.7224960 0.5095933
3_2 3 2 0 1 5.516060 5.314505 5.516060 4.283971 0.5095933 -0.7224960

Distinguishable Dyads - Fitting the Model in BRMS

formula <- bf(
  satisfaction ~ 
    
    # Remove global intercept, introduce male and female intercepts.
    0 + is_male + is_female +
    
    # Actor effect for male and female
    communication_actor_gmc:is_male + 
    communication_actor_gmc:is_female + 
    
    # Partner effect for male and female
    communication_partner_gmc:is_male + 
    communication_partner_gmc:is_female +
    
    # Option 1: Random intercept for male and female (correlated), not identifiable in the cross-sectional case:
    # (0 + is_male + is_female | coupleID)
    
    # Option 2: Too restrictive, basically modelling Exchangeable dyads with perfect correlation
    # (1 | coupleID) 
    
    # Solution: use compound symmetry of the residuals
    cosy(gr = coupleID),
  
  # Two distinct residual variances for males vs females:
  sigma = ~ 0 + is_male + is_female
)

priors <- c(
  # prior(normal(2, 3), class = "Intercept"),
  prior(normal(0, 5), class = "b"),
  prior(student_t(3, 0, 1.5), class = "b", dpar = "sigma"),
  prior(beta(1, 3), class = "cosy")
)

model_dist_apim <- brm(
  formula = formula, 
  data = df_apim,
  family = gaussian(link = identity),
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000, 
  seed = 123,
  file = file.path('brms_cache', 'example1_dist_apim') # Cache the model
)

Distinguishable Dyads - Check Model Convergence and Fit

Check Rhats and Effective Sample Sizes (ESS_tail and ESS_bulk) directly from the brms summary. Additionally you can:

rstan::check_hmc_diagnostics(model_dist_apim$fit)

Divergences:
0 of 4000 iterations ended with a divergence.

Tree depth:
0 of 4000 iterations saturated the maximum tree depth of 10.

Energy:
E-BFMI indicated no pathological behavior.
loo::pareto_k_table(loo(model_dist_apim))

All Pareto k estimates are good (k < 0.7).
plot(model_dist_apim, ask = FALSE)
pp_check(model_dist_apim, 'dens_overlay_grouped', group = 'is_male')
pp_check(model_dist_apim, 'ecdf_overlay_grouped', group = 'is_male')
pp_check(model_dist_apim, type = "loo_pit_overlay")
# Custom function to make DHARMa work with brms (see file 'Functions')
DHARMa.check_brms(model_dist_apim)
Est. 95% CI Rhat Bulk_ESS Tail_ESS
Fixed Effects
is_male 5.63 [5.50, 5.75] 1.002 7151 3296
is_female 4.57 [4.44, 4.70] 1.001 6923 3106
is_male:communication_actor_gmc 1.58 [1.48, 1.67] 1.000 5701 3439
is_female:communication_actor_gmc 1.83 [1.74, 1.91] 1.001 5483 3095
is_male:communication_partner_gmc 0.35 [0.25, 0.44] 1.000 5606 3328
is_female:communication_partner_gmc 0.24 [0.14, 0.33] 1.000 5580 3578
Residual Structure
cosy 0.26 [0.18, 0.34] 1.001 5755 3058
sigma_is_male 0.40 [0.34, 0.46] 1.001 6739 2883
sigma_is_female 0.40 [0.34, 0.46] 1.000 7213 3086

Distinguishable Dyads - (Intensive) Longitudinal APIM

Distinguishable Dyads - L-APIM - Data: Time and Ties

userID coupleID diaryday is_male is_female daily_closeness provided_support_actor provided_support_partner
1063 31 0 1 0 3 0.67 1.33
1063 31 1 1 0 4 1 1.67
1063 31 2 1 0 4 2.67 2
... ... ... ... ... ... ... ...
1064 31 0 0 1 3 1.33 0.67
1064 31 1 0 1 3 1.67 1
1064 31 2 0 1 4 2 2.67

Distinguishable Dyads - L-APIM - Centering

With repeated measures we need to center variables in order to not conflate levels of analysis:

  • Between-Person: The person-mean in relation to the grand mean
  • Within-Person: Daily fluctuations of an individual from their person-mean

While it may seem like we have 3-levels (couple/person/day), by including the means of both partners and their correlations in the model, all information about the couple-level (level 3) is included. This will become clear when comparing the exchangeable APIM to the DIM. We thus use a 2-level model with the dyad as the level of analysis.

Distinguishable Dyads - L-APIM - Centering

df_long_apim <- bmlm::isolate(
  df_long_apim, 
  by = 'userID', # NOT coupleID
  value = c('provided_support_actor', 'provided_support_partner'),
  which = 'both' 
) %>%
  rename(
    provided_support_actor_cwp = provided_support_actor_cw,
    provided_support_partner_cwp = provided_support_partner_cw,
    provided_support_actor_cbp = provided_support_actor_cb,
    provided_support_partner_cbp = provided_support_partner_cb,
  )

Distinguishable Dyads - L-APIM - Model

formula <- bf(
  daily_closeness ~ 
    
    # Intercepts
    0 + is_male + is_female +
    
    # Time Slopes
    is_male:diaryday +
    is_female:diaryday +
    
    # Between-Person APIM
    is_male:provided_support_actor_cbp + 
    is_male:provided_support_partner_cbp +
    
    is_female:provided_support_actor_cbp + 
    is_female:provided_support_partner_cbp +
  
    # Within-Person APIM
    is_male:provided_support_actor_cwp + 
    is_male:provided_support_partner_cwp +
    
    is_female:provided_support_actor_cwp + 
    is_female:provided_support_partner_cwp +
    
    # Accounting for non-independence between partners' means and trajectories 
    # and effect sensitivities via random effects:
    (0 + is_male + is_female +  # Random Intercepts and their correlations
       
       is_male:diaryday + is_female:diaryday +
       
       is_male:provided_support_actor_cwp + 
       is_male:provided_support_partner_cwp +
       is_female:provided_support_actor_cwp + 
       is_female:provided_support_partner_cwp | coupleID ) + 
    
    
    # Accounting for daily non-independence
    
    # Option 1: model correlated residuals on the day with a 
    # compound symmetry structure:
    # cosy(time = diaryday, gr = coupleID),
    
    # Option 2: model daily coupling / common shocks with a random effect:
    (1 | coupleID:diaryday) +
    
    
    # Modelling residuals as autocorrelated within each person (order 1)
    ar(time = diaryday, gr = coupleID:userID, p = 1),
    
  sigma ~ 0 + is_male + is_female   # heteroscedastic residuals
)

priors <- c(
  prior(normal(0, 5), class = "b"),
  prior(student_t(3, 0, 1.5), class = "sd"),
  prior(normal(0.5, 0.3), class = "ar"),
  prior(student_t(3, 0, 1.5), class ="b", dpar = "sigma")
)


model_dist_apim_long <- brm(
  formula = formula, 
  data = df_long_apim,
  family = gaussian(link = identity),
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000, 
  seed = 123,
  file = file.path('brms_cache', 'example1_dist_apim_long') # Cache the model
)

Results

Est. 95% CI Rhat Bulk_ESS Tail_ESS
Fixed Effects
is_male 3.82 [ 3.56, 4.09] 1.004 1285 2238
is_female 3.83 [ 3.57, 4.08] 1.004 1260 2008
is_male:diaryday 0.00 [-0.01, 0.00] 1.000 4295 3592
is_female:diaryday 0.00 [ 0.00, 0.01] 1.002 3938 3105
is_male:provided_support_actor_cbp 0.03 [-0.23, 0.29] 1.000 2608 2848
is_male:provided_support_partner_cbp 0.09 [-0.16, 0.36] 1.002 2307 2518
is_female:provided_support_actor_cbp -0.05 [-0.30, 0.21] 1.001 2626 2875
is_female:provided_support_partner_cbp 0.03 [-0.22, 0.29] 1.001 2679 3109
is_male:provided_support_actor_cwp -0.02 [-0.09, 0.05] 1.001 5130 3362
is_male:provided_support_partner_cwp 0.01 [-0.07, 0.09] 1.000 5303 3387
is_female:provided_support_actor_cwp 0.02 [-0.07, 0.11] 1.001 3999 1896
is_female:provided_support_partner_cwp 0.01 [-0.05, 0.08] 1.002 4152 3168
Random Effects
coupleID: sd(is_male) 0.62 [ 0.45, 0.85] 1.002 2126 2867
coupleID: sd(is_female) 0.60 [ 0.42, 0.81] 1.001 1724 2360
coupleID: sd(is_male:diaryday) 0.01 [ 0.00, 0.02] 1.008 450 726
coupleID: sd(is_female:diaryday) 0.01 [ 0.00, 0.02] 1.018 333 655
coupleID: sd(is_male:provided_support_actor_cwp) 0.06 [ 0.00, 0.19] 1.002 1079 2076
coupleID: sd(is_male:provided_support_partner_cwp) 0.16 [ 0.01, 0.35] 1.005 556 1155
coupleID: sd(is_female:provided_support_actor_cwp) 0.17 [ 0.01, 0.37] 1.015 579 958
coupleID: sd(is_female:provided_support_partner_cwp) 0.26 [ 0.13, 0.44] 1.001 1295 1622
coupleID: cor(is_male,is_female) 0.67 [ 0.34, 0.86] 1.006 1421 2783
coupleID: cor(is_male,is_male:diaryday) -0.13 [-0.62, 0.52] 1.000 2170 2612
coupleID: cor(is_female,is_male:diaryday) 0.10 [-0.47, 0.62] 1.002 2433 2743
coupleID: cor(is_male,is_female:diaryday) 0.03 [-0.54, 0.62] 1.000 3356 2852
coupleID: cor(is_female,is_female:diaryday) -0.08 [-0.59, 0.54] 1.000 2246 1695
coupleID: cor(is_male:diaryday,is_female:diaryday) 0.34 [-0.43, 0.80] 1.005 827 2237
coupleID: cor(is_male,is_male:provided_support_actor_cwp) -0.11 [-0.66, 0.56] 1.001 3669 1130
coupleID: cor(is_female,is_male:provided_support_actor_cwp) -0.22 [-0.71, 0.48] 1.001 3153 1220
coupleID: cor(is_male:diaryday,is_male:provided_support_actor_cwp) -0.02 [-0.63, 0.60] 1.000 3623 3095
coupleID: cor(is_female:diaryday,is_male:provided_support_actor_cwp) 0.01 [-0.62, 0.63] 1.001 3288 3216
coupleID: cor(is_male,is_male:provided_support_partner_cwp) -0.03 [-0.56, 0.52] 1.001 3421 2187
coupleID: cor(is_female,is_male:provided_support_partner_cwp) -0.01 [-0.52, 0.54] 1.000 3786 2964
coupleID: cor(is_male:diaryday,is_male:provided_support_partner_cwp) 0.12 [-0.53, 0.67] 1.002 1579 2797
coupleID: cor(is_female:diaryday,is_male:provided_support_partner_cwp) 0.21 [-0.51, 0.73] 1.003 1271 2489
coupleID: cor(is_male:provided_support_actor_cwp,is_male:provided_support_partner_cwp) 0.05 [-0.60, 0.65] 1.003 1737 2016
coupleID: cor(is_male,is_female:provided_support_actor_cwp) -0.05 [-0.58, 0.49] 1.000 4045 2484
coupleID: cor(is_female,is_female:provided_support_actor_cwp) -0.17 [-0.64, 0.48] 1.002 2733 2200
coupleID: cor(is_male:diaryday,is_female:provided_support_actor_cwp) 0.00 [-0.62, 0.62] 1.001 2263 2945
coupleID: cor(is_female:diaryday,is_female:provided_support_actor_cwp) 0.12 [-0.54, 0.68] 1.002 1860 2809
coupleID: cor(is_male:provided_support_actor_cwp,is_female:provided_support_actor_cwp) 0.14 [-0.57, 0.71] 1.002 1643 1508
coupleID: cor(is_male:provided_support_partner_cwp,is_female:provided_support_actor_cwp) 0.30 [-0.48, 0.80] 1.003 1048 2088
coupleID: cor(is_male,is_female:provided_support_partner_cwp) -0.01 [-0.46, 0.43] 1.001 2655 2871
coupleID: cor(is_female,is_female:provided_support_partner_cwp) -0.27 [-0.64, 0.17] 1.002 2441 2563
coupleID: cor(is_male:diaryday,is_female:provided_support_partner_cwp) -0.13 [-0.66, 0.50] 1.002 1440 2719
coupleID: cor(is_female:diaryday,is_female:provided_support_partner_cwp) 0.08 [-0.55, 0.63] 1.005 1180 2017
coupleID: cor(is_male:provided_support_actor_cwp,is_female:provided_support_partner_cwp) 0.23 [-0.49, 0.77] 1.001 1083 2408
coupleID: cor(is_male:provided_support_partner_cwp,is_female:provided_support_partner_cwp) 0.15 [-0.48, 0.67] 1.000 1793 2065
coupleID: cor(is_female:provided_support_actor_cwp,is_female:provided_support_partner_cwp) 0.22 [-0.42, 0.73] 1.001 1786 2241
coupleID:diaryday: sd(Intercept) 0.56 [ 0.49, 0.62] 1.004 783 2133
Residual Structure
ar[1] 0.09 [-0.01, 0.19] 1.001 1639 2460
sigma_is_male -0.52 [-0.62, -0.43] 1.001 1764 2853
sigma_is_female -0.54 [-0.63, -0.44] 1.001 1638 1933

Exchangeable Dyads

Exchangeable Dyads - Cross-Sectional APIM

Exchangeable APIM

Exchangeable Dyads - Cross-Sectional APIM

Exchangeable APIM

Exchangeable Dyads - Main Assumptions

Partners are exchangeable, i.e., not systematically different.

  • Equal actor effects
  • Equal partner effects
  • Equal means
  • Equal residual variances

But they should still be allowed to vary within each couple, while being correlated.

Exchangeable Dyads - Cross-Sectional APIM: Data

Same cross-sectional data in the same 4-field actor-partner format as before.

print_df(head(df_apim))
userID coupleID gender is_male is_female communication satisfaction communication_actor communication_partner communication_actor_gmc communication_partner_gmc
1_1 1 1 1 0 4.851107 5.214440 4.851107 6.004533 -0.1553604 0.9980657
1_2 1 2 0 1 6.004533 5.812518 6.004533 4.851107 0.9980657 -0.1553604
2_1 2 1 1 0 5.881321 7.467339 5.881321 6.433723 0.8748537 1.4272563
2_2 2 2 0 1 6.433723 7.256714 6.433723 5.881321 1.4272563 0.8748537
3_1 3 1 1 0 4.283971 6.047375 4.283971 5.516060 -0.7224960 0.5095933
3_2 3 2 0 1 5.516060 5.314505 5.516060 4.283971 0.5095933 -0.7224960

Exchangeable Dyads - Cross-Sectional APIM: Model

formula <- bf(
  satisfaction ~ 1 + 
    communication_actor_gmc + communication_partner_gmc + 
    
    # Option 1: A single couple level random intercept. Too restrictive.
    # Would impose perfect correlation between partners in each couple
    # (1 | coupleID)
    
    cosy(gr = coupleID)
  
  # Note: no need to model separate sigmas for each partner.
  # Homogeneous residual variance is estimated: 
  # Implied: sigma = ~ 1
)

priors <- c(
  prior(normal(2, 10), class = "Intercept"),
  prior(normal(0, 5), class = "b"),
  prior(student_t(3, 0, 1.5), class = "sigma"),
  prior(beta(1, 3), class = "cosy")
)

model_ind_apim <- brm(
  formula = formula, 
  data = df_apim,
  family = gaussian(link = identity),
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000, 
  seed = 123,
  file = file.path('brms_cache', 'example1_ind_apim') # Cache the model
)

Exchangeable Dyads - Cross-Sectional APIM: Results

Est. 95% CI Rhat Bulk_ESS Tail_ESS
Fixed Effects
Intercept 5.10 [5.00, 5.21] 1.003 4941 2961
communication_actor_gmc 1.70 [1.64, 1.77] 1.001 4653 3315
communication_partner_gmc 0.29 [0.22, 0.36] 1.000 4476 2728
Residual Structure
cosy 0.13 [0.05, 0.20] 1.002 4388 2183
sigma 1.59 [1.52, 1.66] 1.000 4881 2928

Exchangeable Dyads - Cross-Sectional APIM: Test for Distinguishability

Leave-one-out (loo) cross-validation for model comparison (even if not nested).

a <- loo_compare(
  loo(model_ind_apim), 
  loo(model_dist_apim)
)
print(a)
Model elpd_diff se_diff
model_dist_apim 0.0 0.0
model_ind_apim -98.7 13.0
report::report(a)

The difference in predictive accuracy, as indexed by Expected Log Predictive Density (ELPD-LOO), suggests that ‘model_dist_apim’ is the best model (ELPD = -1965.61), followed by ‘model_ind_apim’ (diff-ELPD = -98.72 ± 13.00, p < .001).
See: documentation of “report”

Exchangeable Dyads - Cross-Sectional DIM

Exchangeable Dyads - Cross-Sectional DIM: Data

Starting from scratch (no 4-field data needed)

userID coupleID communication satisfaction
1_1 1 4.851107 5.214440
1_2 1 6.004533 5.812518
2_1 2 5.881321 7.467339
2_2 2 6.433723 7.256714
3_1 3 4.283971 6.047375
3_2 3 5.516060 5.314505

Exchangeable Dyads - Cross-Sectional DIM: Centering

Decompose communication variance into:

  • Between-couple (cbc): Couple-mean communication skills in relation to other couples
  • Within-couple (cwc): Individuals’ communication skills in relation to their couple-mean

Same assumption about exchangeability as in the APIM

Exchangeable Dyads - Cross-Sectional DIM: Centering

df_dim <- df_apim2 %>%
  group_by(coupleID) %>%
  mutate(
    communication_cm = mean(communication, na.rm = TRUE),
    communication_cwc= communication - communication_cm
  ) %>%
  ungroup() %>%
  mutate(
    communication_cbc= communication_cm - mean(communication_cm, na.rm = TRUE)
  ) %>% 
  select(c('userID', 'coupleID', 'satisfaction', 'communication_cwc', 'communication_cbc'))

print_df(head(df_dim))
userID coupleID satisfaction communication_cwc communication_cbc
1_1 1 5.214440 -0.5767130 0.4213527
1_2 1 5.812518 0.5767130 0.4213527
2_1 2 7.467339 -0.2762013 1.1510550
2_2 2 7.256714 0.2762013 1.1510550
3_1 3 6.047375 -0.6160447 -0.1064514
3_2 3 5.314505 0.6160447 -0.1064514

Exchangeable Dyads - Cross-Sectional DIM: Model

formula <- bf(
  satisfaction ~ 1 + 
    communication_cbc+ communication_cwc + 
    cosy(gr = coupleID) 
)

priors <- c(
  prior(normal(2, 10), class = "Intercept"),
  prior(normal(0, 5), class = "b"),
  prior(student_t(3, 0, 1.5), class = "sigma"),
  prior(beta(1, 3), class = "cosy")
)

model_ind_dim <- brm(
  formula = formula, 
  data = df_dim,
  family = gaussian(link = identity),
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000, 
  seed = 123,
  file = file.path('brms_cache', 'example1_ind_dim') # Cache the model
)

Exchangeable Dyads - Cross-Sectional DIM: Results

Est. 95% CI Rhat Bulk_ESS Tail_ESS
Fixed Effects
Intercept 5.10 [5.00, 5.20] 1.002 4768 2763
communication_cbc 2.00 [1.92, 2.08] 1.005 5242 2876
communication_cwc 1.41 [1.30, 1.52] 1.001 4873 2970
Residual Structure
cosy 0.12 [0.05, 0.20] 0.999 3747 1951
sigma 1.59 [1.52, 1.66] 1.002 4596 2739

Equivalence APIM and DIM

APIM (left) vs. DIM (right)

\[b_{actor\_gmc} + b_{partner\_gmc} = b_{cbc}\] \[b_{actor\_gmc} - b_{partner\_gmc} = b_{cwc}\] (Bolger et al., 2025)

Equivalence APIM and DIM

Results

Top sliders: grand-mean-centered Communication for Actor and Partner (APIM).
Bottom sliders: reparameterized communication — Centered Between-Couple (xcbc) and Centered Within-Couple (xcwc) (DIM).

0.00 0.00 0.00 0.00

Equivalence APIM and DIM

  • If the couple mean goes up by 1 and the within-couple deviation from the couple mean stays fixed, this means that both partners’ predictors must go up by 1 unit. Thus, the effects are linear combinations:

\[b_{cbc} = b_{actor\_gmc} + b_{partner\_gmc}\]

  • Similarly, if a person’s deviation from their couple mean is one unit higher and the couple mean stays fixed, this means their partners’ deviation must be one unit lower. Thus, the effects are a linear combination:

\[b_{cwc} = b_{actor\_gmc} - b_{partner\_gmc}\]

Equivalence APIM and DIM

  • Conversely, to obtain the actor and partner effects given the DIM estimates:

\[b_{actor\_gmc} = \frac{b_{cbc} + b_{cwc}}{2}\]

\[b_{partner\_gmc} = \frac{b_{cbc} - b_{cwc}}{2}\]

Equivalence APIM and DIM

Example: using APIM coefficients to obtain DIM coefficients and computing Credible Intervals of DIM coefficients.

a <- hypothesis(
  model_ind_apim, 
  "communication_actor_gmc + communication_partner_gmc = 0"
)

a$hypothesis[,c(2,3,4,5)]
  Estimate Est.Error CI.Lower CI.Upper
1 1.995968 0.0411011 1.915506 2.076755

APIM (left) vs. DIM (right)

Equivalence APIM and DIM: Takeaway

  • APIM and DIM are reparametrizations of the same model
  • APIM: intuitive actor/partner framing
  • DIM: clean separation of between vs within components
  • Estimating one model allows for directly obtaining estimates of the other
  • Same random-effects structure at dyad level and same assumptions

Side Note:

In the distinguishable case, the Dyadic Score Model (DSM) (e.g., Stadler et al., 2023) is equivalent to the APIM! (Bolger et al., 2025; Iida et al., 2018)

Exchangeable Dyads - Longitudinal APIM/DIM: Data

df_long_apim2 <- df_long_apim %>%
  select('coupleID', 'userID', 'diaryday', 'daily_closeness', 
         'provided_support_actor_cwp', 'provided_support_partner_cwp', 
         'provided_support_actor_cbp', 'provided_support_partner_cbp')

print_df(print_couple_preview(df_long_apim2, couple_id = "31"))  
coupleID userID diaryday daily_closeness provided_support_actor_cwp provided_support_partner_cwp provided_support_actor_cbp provided_support_partner_cbp
31 1063 0 3 -0.32 0.73 0.18 -0.19
31 1063 1 4 0.02 1.06 0.18 -0.19
31 1063 2 4 1.68 1.39 0.18 -0.19
... ... ... ... ... ... ... ...
31 1064 0 3 0.73 -0.32 -0.19 0.18
31 1064 1 3 1.06 0.02 -0.19 0.18
31 1064 2 4 1.39 1.68 -0.19 0.18

Exchangeable Dyads - Longitudinal DIM/APIM

  • We need an APIM or DIM on both the between person level and the within-person level.
  • Due to equivalence, we could have a DIM on one level and an APIM on the other.

If we want a within-person DIM:

df_long_dim <- df_long_apim2 %>%
  mutate(
    
    # Between Person Level DIM (centering between couples and within couple between person)
    provided_support_cbc = (provided_support_actor_cbp + provided_support_partner_cbp) / 2,
    provided_support_cwcbp = (provided_support_actor_cbp - provided_support_partner_cbp) / 2,
    
    # Within Person Level DIM (couple mean and deviation on each day)
    provided_support_cwp_mean = (provided_support_actor_cwp + provided_support_partner_cwp) / 2,
    provided_support_cwp_halfdiff = (provided_support_actor_cwp - provided_support_partner_cwp) / 2
  ) 

Exchangeable Dyads - Longitudinal DIM/APIM: Preparing the sum and difference approach

We need to randomly assign -1 and +1 for each Partner within each couple. This will be needed for contrast coding the intercept.

For example:

df_long_dim <- df_long_dim %>%
  group_by(coupleID) %>%
  mutate(
    base = ifelse(userID == min(userID), 1, -1),
    flip = 1 - 2 * rbinom(1, 1, 0.5),   # yields +1 or -1 once per couple
    Idiff = base * flip
  ) %>%
  ungroup() %>%
  relocate(Idiff, .after = userID) %>%
  select(-base, -flip)


print_df(print_couple_preview(df_long_dim, couple_id = "31"))  

Exchangeable Dyads - Longitudinal DIM/APIM: Preparing the Sum and Difference Approach

coupleID userID Idiff diaryday daily_closeness provided_support_actor_cwp provided_support_partner_cwp provided_support_actor_cbp provided_support_partner_cbp provided_support_cbc provided_support_cwcbp provided_support_cwp_mean provided_support_cwp_halfdiff
31 1063 -1 0 3 -0.32 0.73 0.18 -0.19 -0.01 0.19 0.2 -0.52
31 1063 -1 1 4 0.02 1.06 0.18 -0.19 -0.01 0.19 0.54 -0.52
31 1063 -1 2 4 1.68 1.39 0.18 -0.19 -0.01 0.19 1.54 0.14
... ... ... ... ... ... ... ... ... ... ... ... ...
31 1064 1 0 3 0.73 -0.32 -0.19 0.18 -0.01 -0.19 0.2 0.52
31 1064 1 1 3 1.06 0.02 -0.19 0.18 -0.01 -0.19 0.54 0.52
31 1064 1 2 4 1.39 1.68 -0.19 0.18 -0.01 -0.19 1.54 -0.14

Exchangeable Dyads - Longitudinal APIM/DIM: Model

formula <- bf(
  daily_closeness ~ 1 + 
    
    diaryday +
  
    # Within-person APIM
    provided_support_actor_cwp + 
    provided_support_partner_cwp +
    # Equivalent to within-person DIM:
    # provided_support_cwp_mean + 
    # provided_support_cwp_halfdiff +
  
    # Between-person APIM
    provided_support_actor_cbp +
    provided_support_partner_cbp +
    # Equivalent to between-person DIM
    # provided_support_cwcbp+
    # provided_support_cbc +
    
    # Dyad-Level intercept and slopes for time-varying predictors
    (1 + diaryday + provided_support_actor_cwp + provided_support_partner_cwp | coupleID) +
    # Both partners' deviations from these dyad-level means and slopes
    (0 + Idiff + 
         I(Idiff * diaryday) + 
         I(Idiff * provided_support_actor_cwp) + 
         I(Idiff * provided_support_partner_cwp) | coupleID) +
    
    # Accounting for same-day shocks/coupling
    (1 | coupleID:diaryday) +
    # alternatively and equivalently if no autocorrelation needed:
    # cosy(time = diaryday, gr = coupleID) 
    
    # Autocorrelated residuals
    ar(time = diaryday, gr = coupleID:userID, p = 1)
  
  # Again, no need to model heterogeneous residual variances (sigma)
  # Implied: sigma = ~ 1
)

priors <- c(
  prior(normal(2, 10), class = "Intercept"),
  prior(normal(0, 5), class = "b"),
  prior(student_t(3, 0, 1.5), class = "sd"),
  prior(lkj(2), class = "cor"), # correlation prior for random effect matrix
  prior(normal(0.5, 0.3), class = "ar"),
  prior(student_t(3, 0, 1.5), class = "sigma")
)

model_apim_ind_long <- brm(
  formula = formula, 
  data = df_long_dim,
  family = gaussian(link = identity),
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000, 
  seed = 123,
  file = file.path('brms_cache', 'model_apim_ind_long_both')
)

Exchangeable Dyads - Longitudinal APIM/DIM: Explanation Idiff

For appropriate random effects, we can use the Sum and Difference approach (del Rosario & West, 2025; Kenny & Ackerman, 2023):

  • Randomly give one partner a constant of -1 and the other a constant of 1
  • The couple level intercept represents the mean (or sum) of both partners’ intercepts. (1 | coupleID)
  • A column of 1s and -1s represent deviations (difference) of each partner from the couple intercept with each partner contributing equally in one direction (+1 and -1)

Exchangeability is retained: if partners are flipped, the result is the same.

(del Rosario & West (2025) provide practical guidance on how to reduce the random effects structure in case of non-convergence.)

Exchangeable Dyads - Longitudinal APIM/DIM: Results

Comparing APIM (left) and DIM (right) output

APIM (left) vs. DIM (right) APIM

Equivalence holds on both levels

Exchangeable Dyads - Longitudinal APIM/DIM: Results

APIM Results in Detail:

Est. 95% CI Rhat Bulk_ESS Tail_ESS
Fixed Effects
Intercept 3.85 [ 3.58, 4.11] 1.001 987 1546
diaryday 0.00 [ 0.00, 0.00] 1.000 3834 3453
provided_support_actor_cwp -0.01 [-0.06, 0.04] 1.001 3425 3111
provided_support_partner_cwp 0.01 [-0.04, 0.07] 1.001 3933 2824
provided_support_actor_cbp -0.01 [-0.15, 0.13] 1.001 2643 2960
provided_support_partner_cbp 0.04 [-0.10, 0.18] 1.001 2721 2721
Random Effects
coupleID: sd(Intercept) 0.66 [ 0.47, 0.92] 1.001 1182 1935
coupleID: sd(diaryday) 0.01 [ 0.01, 0.02] 1.007 652 677
coupleID: sd(provided_support_actor_cwp) 0.14 [ 0.03, 0.25] 1.002 667 577
coupleID: sd(provided_support_partner_cwp) 0.23 [ 0.13, 0.35] 1.001 983 1296
coupleID: sd(Idiff) 0.17 [ 0.10, 0.27] 1.003 1413 1345
coupleID: sd(IIdiffMUdiaryday) 0.00 [ 0.00, 0.01] 1.002 877 1867
coupleID: sd(IIdiffMUprovided_support_actor_cwp) 0.05 [ 0.00, 0.15] 1.005 1069 1736
coupleID: sd(IIdiffMUprovided_support_partner_cwp) 0.11 [ 0.01, 0.22] 1.004 926 1130
coupleID: cor(Intercept,diaryday) -0.33 [-0.69, 0.23] 1.002 1259 2024
coupleID: cor(Intercept,provided_support_actor_cwp) -0.37 [-0.78, 0.25] 1.003 1205 1582
coupleID: cor(diaryday,provided_support_actor_cwp) 0.28 [-0.35, 0.75] 1.007 1302 2040
coupleID: cor(Intercept,provided_support_partner_cwp) -0.29 [-0.68, 0.20] 1.003 1498 2460
coupleID: cor(diaryday,provided_support_partner_cwp) 0.32 [-0.25, 0.74] 1.001 1331 1805
coupleID: cor(provided_support_actor_cwp,provided_support_partner_cwp) 0.63 [-0.05, 0.92] 1.003 867 813
coupleID: cor(Idiff,IIdiffMUdiaryday) -0.17 [-0.78, 0.61] 1.001 3228 3067
coupleID: cor(Idiff,IIdiffMUprovided_support_actor_cwp) 0.09 [-0.62, 0.73] 1.002 3147 2747
coupleID: cor(IIdiffMUdiaryday,IIdiffMUprovided_support_actor_cwp) 0.04 [-0.68, 0.73] 1.001 2779 3004
coupleID: cor(Idiff,IIdiffMUprovided_support_partner_cwp) -0.27 [-0.76, 0.40] 1.001 2244 2573
coupleID: cor(IIdiffMUdiaryday,IIdiffMUprovided_support_partner_cwp) 0.03 [-0.69, 0.71] 1.003 1196 2093
coupleID: cor(IIdiffMUprovided_support_actor_cwp,IIdiffMUprovided_support_partner_cwp) -0.14 [-0.79, 0.61] 1.002 1413 2300
coupleID:diaryday: sd(Intercept) 0.53 [ 0.47, 0.60] 1.001 776 1769
Residual Structure
ar[1] 0.10 [ 0.00, 0.19] 1.000 1566 2402
sigma 0.59 [ 0.56, 0.63] 1.003 1241 1745

Extract Full APIM Random Effects Variance-Covariance Matrix

We can rotate the random effects structure back to obtain the full APIM Random effects variance-covariance matrix. Just as we would in a SEM model with equality constraints!

# Custom function 
rotate_apim_covariance <- function(
    fit,
    gr = "coupleID",
    Idiff = "Idiff"
) {
  random_effects <- fit$ranef
  random_effects <- random_effects[random_effects$group == gr,]

  SUM <- random_effects$coef[!grepl(Idiff, random_effects$coef)]
  DIFF <- random_effects$coef[grepl(Idiff, random_effects$coef)]
  
  varcor <- VarCorr(fit)[[gr]]
  
  within_person_covariance_matrix <- varcor$cov[SUM,'Estimate',SUM] + varcor$cov[DIFF,'Estimate',DIFF]
  cross_person_covariance_matrix <- varcor$cov[SUM,'Estimate',SUM] - varcor$cov[DIFF,'Estimate',DIFF]
  
  p <- nrow(within_person_covariance_matrix)
  
  Full <- rbind(
    cbind(within_person_covariance_matrix,  cross_person_covariance_matrix),
    cbind(cross_person_covariance_matrix, within_person_covariance_matrix)
  )
  
  # nice labels
  base <- SUM
  rn <- c(paste0("PartnerA_", base), paste0("PartnerB_", base))
  colnames(Full) <- rn
  rownames(Full) <- rn
  
  return(list(within_person_covariance_matrix=within_person_covariance_matrix, cross_person_covariance_matrix=cross_person_covariance_matrix, full_covariance_matrix=Full))
}

a <- rotate_apim_covariance(model_apim_ind_long)

Extract Full APIM Random Effects Variance-Covariance Matrix: Within-Person Matrix

print_df(a$within_person_covariance_matrix)
Intercept diaryday provided_support_actor_cwp provided_support_partner_cwp
Intercept 0.4999548 -0.0028921 -0.0354603 -0.0505549
diaryday -0.0028921 0.0001557 0.0004753 0.0008553
provided_support_actor_cwp -0.0354603 0.0004753 0.0281550 0.0195700
provided_support_partner_cwp -0.0505549 0.0008553 0.0195700 0.0709432

Extract Full APIM Random Effects Variance-Covariance Matrix: Cross-Person Matrix

print_df(a$cross_person_covariance_matrix)
Intercept diaryday provided_support_actor_cwp provided_support_partner_cwp
Intercept 0.4334127 -0.0027084 -0.0373154 -0.0394890
diaryday -0.0027084 0.0001458 0.0004602 0.0008405
provided_support_actor_cwp -0.0373154 0.0004602 0.0186953 0.0226297
provided_support_partner_cwp -0.0394890 0.0008405 0.0226297 0.0402293

Extract Full APIM Random Effects Variance-Covariance Matrix: Full Matrix

print_df(a$full_covariance_matrix)
PartnerA_Intercept PartnerA_diaryday PartnerA_provided_support_actor_cwp PartnerA_provided_support_partner_cwp PartnerB_Intercept PartnerB_diaryday PartnerB_provided_support_actor_cwp PartnerB_provided_support_partner_cwp
PartnerA_Intercept 0.4999548 -0.0028921 -0.0354603 -0.0505549 0.4334127 -0.0027084 -0.0373154 -0.0394890
PartnerA_diaryday -0.0028921 0.0001557 0.0004753 0.0008553 -0.0027084 0.0001458 0.0004602 0.0008405
PartnerA_provided_support_actor_cwp -0.0354603 0.0004753 0.0281550 0.0195700 -0.0373154 0.0004602 0.0186953 0.0226297
PartnerA_provided_support_partner_cwp -0.0505549 0.0008553 0.0195700 0.0709432 -0.0394890 0.0008405 0.0226297 0.0402293
PartnerB_Intercept 0.4334127 -0.0027084 -0.0373154 -0.0394890 0.4999548 -0.0028921 -0.0354603 -0.0505549
PartnerB_diaryday -0.0027084 0.0001458 0.0004602 0.0008405 -0.0028921 0.0001557 0.0004753 0.0008553
PartnerB_provided_support_actor_cwp -0.0373154 0.0004602 0.0186953 0.0226297 -0.0354603 0.0004753 0.0281550 0.0195700
PartnerB_provided_support_partner_cwp -0.0394890 0.0008405 0.0226297 0.0402293 -0.0505549 0.0008553 0.0195700 0.0709432

References

Bolger, N., Laurenceau, J.-P., & DiGiovanni, A. (2025). Unified analysis model for indistinguishable and distinguishable dyads. Innovations in Interpersonal Relationships and Health Research: Advancing the Integration of Interdisciplinary Approaches to Dyadic Behavior Change. https://doi.org/10.17605/OSF.IO/WYDCJ
del Rosario, K. S., & West, T. V. (2025). A Practical Guide to Specifying Random Effects in Longitudinal Dyadic Multilevel Modeling. Advances in Methods and Practices in Psychological Science, 8(3), 25152459251351286. https://doi.org/10.1177/25152459251351286
Iida, M., Seidman, G., & Shrout, P. E. (2018). Models of interdependent individuals versus dyadic processes in relationship research. Journal of Social and Personal Relationships, 35(1), 59–88. https://doi.org/10.1177/0265407517725407
Kenny, D. A., & Ackerman, R. A. (2023). Estimation of random effects with over-time dyadic data using multilevel modeling: The sum and difference method. OSF Preprints. https://doi.org/10.31219/osf.io/fju72
Stadler, G., Scholz, U., Bolger, N., Shrout, P. E., Knoll, N., & Lüscher, J. (2023). How is companionship related to romantic partners’ affect, relationship satisfaction, and health behavior? Using a longitudinal dyadic score model to understand daily and couple-level effects of a dyadic predictor. Applied Psychology: Health and Well-Being, 15(4), 1530–1554. https://doi.org/10.1111/aphw.12450